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A number of theoretically well-motivated additions to the standard cosmological model predict 
weak signatures in the form of spatially localized sources embedded in the cosmic microwave back- 
ground (CMB) fluctuations. We present a hierarchical Bayesian statistical formalism and a complete 
data analysis pipeline for testing such scenarios. We derive an accurate approximation to the full 
posterior probability distribution over the parameters defining any theory that predicts sources em- 
bedded in the CMB, and perform an extensive set of tests in order to establish its validity. The 
approximation is implemented using a modular algorithm, designed to avoid a posteriori selec- 
tion effects, which combines a candidate-detection stage with a full Bayesian model-selection and 
parameter-estimation analysis. We apply this pipeline to theories that predict cosmic textures and 
bubble collisions, extending previous analyses by using: (1) adaptive-resolution techniques, allowing 
us to probe features of arbitrary size, and (2) optimal filters, which provide the best possible sensi- 
tivity for detecting candidate signatures. We conclude that the WMAP 7-year data do not favor the 
addition of either cosmic textures or bubble collisions to ACDM, and place robust constraints on the 
predicted number of such sources. The expected numbers of bubble collisions and cosmic textures 
on the CMB sky are constrained to be fewer than 4.0 and 5.2 at 95% confidence, respectively. 
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I. INTRODUCTION 



The cosmic microwave background (CMB) radiation provides our best picture of the primordial universe, and 
therefore the best set of observations available to confront theories of the early universe with data. The angular 
power spectrum of the CMB, together with complementary datasets (e.g., from large-scale structure and supernova 
surveys), has established the standard model of cosmology, a spatially flat universe dominated by cold dark matter 
(ACDM). However, many theories of high-energy physics predict that there should be deviations from the isotropic 
and purely Gaussian density fluctuations predicted by ACDM. In this paper, we are concerned with the question of 
how to optimally test theories that predict spatially-localized sources embedded in the CMB. We present a statistical 
formalism and a set of approximations that are implemented in a full analysis pipeline to construct the posterior 
probability distribution over the parameters describing a class of theories. We implement a two-step algorithm in 
which we first locate the most promising candidate signatures, and then use the information about the number, 
location, and properties of the candidate sources to construct an approximation to the full posterior probability 
distribution. 

To illustrate the application of this pipeline, we focus on two signatures that are predicted by theories with spon- 
taneous symmetry breaking giving rise to phase transitions in the early universe: cosmic textures and cosmic bubble 
collisions. Cosmic textures are a type of topological defect produced when a non-Abelian global symmetry is bro- 
ken pp. Textures are not stable, but instead undergo collapse as they come within the expanding cosmological 
horizon, eventually unwinding into scalar radiation [I}{5]. CMB photons passing through a collapsing texture will be 
red-shifted, while those passing through an exploding texture will be blue-shifted, giving rise to a set of features in 
the CMB [5] . Cosmic bubble collisions are predicted by theories of eternal inflation, where our observable universe 
is postulated to be embedded inside one bubble among many, formed during a first-order phase transition out of an 
inflating false vacuum (for a review of eternal inflation see, e.g., Refs. [3 [7]). The density perturbations induced by 
collisions between our bubble and others can lead to localized features in the CMB, providing an observable signature 
of the dynamics of eternal inflation [8] . 

In previous work, we presented the first constraints on theories giving rise to cosmic bubble collisions [HI 110] and 
the first full-sky constraints on cosmic textures [11] . The present paper focuses on: 

• Generalizing the statistical formalism and approximation scheme used in Refs. [SHU]; 



Implementing an adaptive-resolution analysis, allowing us to overcome the limitations in Refs. [TO] on the 
size of candidate bubble collisions; 

Including and refining the candidate detection scheme using optimal filters presented in Ref. |12j : 

Performing a complete suite of tests of the formalism, approximations, and analysis pipeline; 

Performing a new analysis of the posterior probability distribution for bubble collisions and cosmic textures that 
includes new candidates from the optimal filtering step in combination with the upgraded adaptive-resolution 
analysis pipeline. 



The paper is structured as follows. In Sec. [TTJ we outline the formalism and approximations we use. In Sec. |III[ 
we describe the theoretical predictions for cosmic textures and bubble collisions. The algorithm used to calculate 
the approximated posterior is described in Sees. |IV|V] and tested in Sees. |VI| and |VII[ A null test of the pipeline 



is carried out in Sec. VIII before the pipeline is applied to CMB data from the Wilkinson Microwave Anisotropy 



Probe (WMAP) [13] in Sec. IX The results of this analysis are compared with previous analyses in Sec. [XJ and our 
conclusions are summarised in Sec. IXII 



II. HIERARCHICAL BAYESIAN SOURCE DETECTION FORMALISM 



A. The theory 



The observed fluctuations in the CMB can be modeled as the realization of a random field on the sphere, which, 
under the assumption of isotropy and Gaussianity, is completely characterized by its angular power spectrum. A 
number of extensions of this model predict various populations of distinct sources embedded in the background 
random field. This includes astrophysical sources such as clusters of galaxies (which affect the CMB through the 
Sunyaev-Zel'dovich effect [14] ), and primordial sources such as cosmic textures and cosmic bubble collisions. We 
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restrict our attention to cases where the temperature anisotropies can be described as 



AT 

— (9, 0) = 8(9, 0) + n(9, <f>) + J2 <>>)■ 



(1) 



Here, 8(9, <f) is a realization of the background random field, n(9, 4>) describes the instrumental noise as well as residual 
foregrounds, and ti(9, 4>) are templates for the temperature anisotropies laid down by each of N s distinct sources. All 
terms other than the instrumental noise are assumed to include the effects of a finite instrumental beam. Such a 
theory can be described by 

• Model parameters: This includes the parameters describing the background random field and the source 
templates)]] The background random field is described by the parameters of ACDM, which we denote by the 
vector mACDM- These parameters include: the fraction of energy density in baryons, f2f,/i 2 ; cold dark matter, 
Ocdm^ 2 ; and dark energy, fl^; the scalar spectral index, n s ; the primordial scalar amplitude, A s ; and the optical 
depth to rcionization, r. Modeling the instrument gives a characterization of the expected noise properties. No 
model of the Galactic foreground residuals is available for the dataset considered in this analysis, and we therefore 
resort to null tests of simulations including foreground residuals in order to determine their effects (although a 
model of the foreground residuals could, in principle, be included in the formalism). 

It is convenient to treat the extension hypothesis as a hierarchical Bayesian model (e.g. Ref. [15] ) in which the 
population level parameters are considered separately from the lower-level parameters describing the individual 
sources. The parameters describing the templates are hence divided into two categories: global parameters, m , 
which describe the source population as a whole; and local parameters, m^, characterizing individual sources. 
Any model will possess at least one global parameter N s , the expected total number of detectable sources - 
in addition to any properties common to all templates. Further, any model will possess at least one set of local 
parameters: {9i,(j>i}, the central position of the i th template. Other properties that can differ from template 
to template (e.g., size) are also classified as local parameters. Global template parameters, in addition to the 
parameters of ACDM, can be thought of as labeling different theories, characterizing the background cosmology 
and the type of source. Local parameters characterize the properties of sources in the context of a specific 
theory. 

• Theoretical priors: An important component of the theory is the prior probability distribution over the 
model parameters. In principle, a complete theory of cosmology would provide an explanation for the observed 
properties of the population of sources and the background random field. In general such a full theory is 
not available. To make progress, we will assume that there are no correlations between the properties of the 
background field and the properties of the sources, rendering the prior separable. In the context of a specific 
theory of sources, the prior over local parameters can be fully determined. The priors over the parameters of 
ACDM and the global template parameters are somewhat less certain in the absence of a theoretical construction 
in which different models can be compared^] A reasonable assumption is therefore to use an uninformative prior, 
which assigns equal probability to all possibilities. The use of uninformative priors requires care to be taken 
when defining the prior range, and will be discussed in later sections. 

• Model statistics: Given a set of model parameters, it is necessary to understand how particular realizations of 
the temperature anisotropies are determined. For the sources, this is most efficiently encoded in the theoretical 
prior over local model parameters of the templates. For the background random field and instrumental noise, this 
is most efficiently encoded in the two-point correlation function. Under ACDM, for perfect data, the correlation 
in the temperature between two positions on the sphere is given by 



where 9ij is the angular distance between two points on the sphere labeled by i and j, and Cg (hiacdm) is the 
angular power spectrum, which is dependent on the choice of parameters miacdm- The characterization of the 
instrumental noise and beam depends on the experiment in question, and will be described in greater detail 
below for the WMAP experiment. 



1 The formalism could also be extended to allow for marginalization over any imperfectly known experimental parameters. For simplicity, 
we assume that the parameters of the WMAP experiment are perfectly known. 

2 Of course, the best example is the eternally inflating multiverse, in which regions with diverse physical properties are sampled. Defining 
the theoretical prior in this case is difficult due to the infinite number of regions that must be compared; this is the "Measure Problem" 
of eternal inflation (see Ref. |16| for a recent review) . 




(2) 



t 
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B. The full posterior 



Having fully specified the theory, we can now ask how to test it. Our goal is to construct the posterior probability 
distribution over the global source parameters, given a dataset d consisting of pixclizcd temperature measurements 
covering a solid angle £l bs — 47r/ s k y of the sky (and, optionally, any statistics derived from them). Bayes' theorem 
gives the posterior as 

Prfm Id f \- p r( m o) p r(d|m ,/ sky ) 

Pr(m |d,/ Sky ) - , (3) 

where Pr(mo) is the prior distribution on the global parameters, mo, Pr(d|mo, /sky) is the likelihood of getting the 
observed data, and Pr(d|/ sky ) ensures that the posterior is normalized. The posterior can also be used to derive 
summary statistics, such as confidence intervals on the global model parameters. The quantity N s is always included 
in the set of global model parameters mo, and ACDM is specified by N s = 0. Therefore, we can also perform model 
selection by comparing the posterior probability of a model for which N s = and one which admits N s > 0. 
The model likelihood Pr(d|mo, /sky) is obtained by marginalizing over: 

1. The parameters of ACDM; 

2. The actual number of sources present on the observable sky (given the expected number of sources, N s , the 
actual number is a realization of a Poisson-like process of mean / s k y A^ s ); 

3. The local template parameters. 
Unpacking the model likelihood, we therefore have: 

™ (f, W )N S -U y N s 

Pr(d|m , /sky ) = £ f (4) 

x J dm ACDM Pr(mACDM) / dmi . . . dmj\r s Pr(mi, . . . mjv 8 )Pr(d|iVs, / sky , hiacdm, mo, mi, . . . m Ns ), 

where Pr(rriACDM) is the prior over the parameters of ACDM, Pr(m 1; . . . m^rj is the prior over the local model 
parameters for N s sources, and Pr(d|iV s , / s k yj ^acdMi mo, mi, ■ ■ ■ mjv a ) is the likelihood. For measurements of the 
CMB temperature anisotropics under the assumption of ACDM as the background random field, the likelihood is 
written as 

Pr(d|iVs,/sk y ,m ACD M,mo ! m 1 ,...m Ar J = ^y^rj^, e -( d ~S& t^C^d-^ t(m,)) T /2 ; (5) 

where C is the pixel-pixel covariance matrix including a ACDM CMB signal and instrumental noise. 

One must evaluate Eq. [5] in order to construct the full posterior, Eq. [3] Evaluating this expression is impossible, 
as it requires marginalizing over a formally infinite dimensional parameter space. Even if the parameter space were 
finite, the enormous size of modern CMB datasets, such as those produced by the WMAP or Planck [T7] experiments, 
makes inverting the (non-diagonal) covariance matrix prohibitively expensive. Nevertheless, it is possible to apply a 
controlled and testable sequence of approximations to estimate the full posterior, as we now describe. 



C. Approximation to the posterior 

In order to proceed, we must construct a suitable approximation to the model likelihood. Let us first deal with 
the cosmological parameters niACDM- We assume that our inferences on do not vary significantly for the range of 
ACDM parameters allowed by current cosmological data. Under this assumption, we can proceed as if Pr(mACDM) = 
<5(rriACDM — Aacdm)) where hiacdm is the best-fit concordance cosmological model. Performing the integral over 
mACDM we obtain the approximation to the model likelihood 

Pr(d|m , / sky ) c ^ Pr(d|Ag, (6) 

7V S — s ' 

where we define the helpful short-hand 



Pr(d|iVs) = / dmi . . . dmjv s Pr(mi, . . . m A r s )Pr(d|A r s, / s k y , m A cDM, m , mi, . . . mjvj. 



(7) 



likelihood zero 




FIG. 1. A schematic depicting the approximation scheme we employ. By locating a set of candidate sources, it is possible 
to determine in which regions of parameter space the likelihood function is appreciably different from zero. Eq. [7] can be 
approximated by integrating over only those regions where the likelihood is large. This in addition collapses the sum in Eq. [6] 
to a finite number of terms. Finally, we neglect correlations of the random Gaussian background CMB between pixels inside 
and outside each blob. 



In a similar spirit, if we knew something about the dependence of the likelihood on the local model parameters, 
it would be possible to approximate the remaining integrals. This is depicted schematically in Fig. [I] To see how 
this works in detail, let us begin with a particular example. Imagine that there is a region of the sky which has 
been judged by some independent method to be a good candidate source. We can segment the data into a "blob" , 
containing the candidate source, and the rest of the sky. The details of the size and shape of the blob are treated in 
abstract here, and will depend on the particular theory of the sources being tested. We can now evaluate the sum 
over N s in Eq. [6] term by term. The likelihood in the first term, for iV s = 0, is simply given by 

Pr(d|iV s = 0) = , _^ /2|c| e- dC ~ ldT/2 , (8) 



(2tt 



which is the likelihood for the null hypothesis, i.e., no sources. Here, and in what follows, C is evaluated at the best 
fit values SiacdM' Moving on to the N a = 1 term, we must evaluate the integral over mi. Recall that the local model 
parameters always include the location at which the template is centered and, if relevant, its size. We can therefore 
separate the integral over mi into the region inside the blob containing our candidate source, which we will refer to 
as region b, and the rest of the sky, which we will refer to as region b: 



Pr(d|iV s = 1) =^dm 1 Pr(m 1 ) ^ r^"^ 



/ dmiPr(mi) (2^ 



-(d-t)C 1 (d-t) T /2 



( 27r )iW2| C | 



(9) 



If there are no sources in the region 6, then we can approximate the likelihood by integrating over region b alone. 
The accuracy of this approximation depends on our ability to locate the candidate source; however, it will always 
provide a lower (i.e., conservative) bound on the likelihood (since we are integrating a positive-definite function). We 
therefore have 



Pr(d|JV s = 1) ~ / d mi Pr( mi ) . 1 ,„,_, e-( d - t > C ~ 1 ( d - t > T / 2 



</ 2 |C| 



(10) 



ib (^y 

While we have reduced the parameter space over which we must integrate, this expression is still numerically intractable 
for large datasets, since we must invert the large, generally non-diagonal covariance matrix. To make progress, we 
must make a few further approximations. Expanding the exponential in the likelihood, we have 

1 -(d-t)C- 1 (d-t) T /2 = 1 r -(d l -U)C'- j 1 (d i -t i )/2 v p -d a C'Jd f3 /2 x e _(d i -t l )C l - 1 d Qi / ir . 



( 27r )AW2|C| 



(2tt 



</3|C 



e -(d l -U)C- J 1 (d J -t J )/2 x e -d a C-£d e /2 
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where the indices i and j refer to pixels in region b while the indices a and /3 correspond to pixels in region b. We 
have used the fact that the template vanishes in region b. Re-arranging, we obtain 



1 P -(d-t)C (d-t) T /2 = e ^ ( 1 p -dC d T /2 



(d-t)C < "' ) 1 (d-t) T /2 / i 



e -dC <66) "dT/2 
-(d-t)C <i ' i ' ) ~ 1 (d-t) T /2 



In these expressions, is the covariance matrix constructed using only the data in region b. We have made two 
approximations in deriving this final expression. First, we have neglected correlations between the template and the 
data in region b, which is equivalent to assuming 

e^UC^dl „ L (15) 

This is justified in the limit where the inverse covariance falls off sufficiently fast with angular distance. Our second 
approximation was to assume that we can make the replacement 

(di - t^Cr^dj - tt) ->• (d- t)C^\d - t) T /2, (16) 

which is justified to the extent that the subgroup of the inverse of the full covariance matrix corresponding to pixels 
in region b can be approximated as the inverse of a covariance matrix defined only using pixels in region b. For a 
diagonal covariance, this is exact. For the non-diagonal covariance matrix of ACDM it is only approximate. We 
comment on the validity of these approximations under ACDM in Section |VII| 
Finally, performing the integral over mi in Eq. 10 we obtain 

Pr(d| AT S = 1) ~ Pr(d|iV s = 0) p 6 (m ), (17) 

where the patch-based evidence ratio pb is given by 

f J dm 1 Pr(m 1 ) e -( d - t ) C(!>b) ~ 1 ( d - t ) T / 2 

This is a measure of how much better the theory with a template at fixed mo fits the patch of data than the theory 
with only the random field. 

Now, we evaluate the two-source term. Again, we approximate the full integral as the integral over region b alone. 
This yields 

Pr(d|Af s = 2)= / /dm 1 dm 2 Pr(m 1 )Pr(m 2 ) / * m -(d-t(m 1 )-t(m 2 ))C- 1 (d-t(m 1 )-t(m 2 )) T /2^ ^ ig j 
JbJb (27t) jv p-/^C| 

Making the same approximation about the covariance matrix as above, we have 

-(66)— 1 



J b J b dmidm 2 Pr(mi)Pr(m 2 )e-( d - t ( mi )- t ( m2 )' C (d-t( mi )-t(m 2 )) T /2 



Pr(d|jV s = 2)^Pr(d|jV s = 0) JfcJb ^~„--v~*,- (2Q) 

g— dC d T /2 

If there is in fact only a single source in region 6, then the addition of another template will not increase the likelihood. 
In effect, we would be trying to fit a single feature with a template possessing twice the number of parameters (mi and 
m 2 ). This will introduce a natural "Occam factor" that favors the simpler model (i.e., the model with one template). 
As a concrete example, assume that the source location is the only local model parameter. If the source can be located 
anywhere on the sky with equal probability, the theory prior is simply 

Pr(m,) = -L. (21) 

In the case where the likelihood function is roughly equal for all positions inside b, with the solid angle contained in b 
given by fif,, and there is no improvement from adding a second template, the relative size of the N s = 1 and N s = 2 
terms can be estimated as 

Pr(d|iVs = 1) „ 4^ 

Pr(d|A^ s = 2) JV 1 ' 
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Assuming the blob does not cover the entire sky, this is always larger than one. Subsequent terms in the N s expansion 
will be penalized by higher powers of this ratio. While this is a toy model, this property is expected to hold generally. 
For a single blob, we can therefore approximate the full-sky posterior Eq. [3] as 

p , , , . s Pr(m ) Pr(d| JV B = 0) (l + (/ sky Ag p b (m )) 

Pr(m |d, /sky ) _ l>mJ - ) , (23) 

where 

Pr(d|/ sky ) = J dm Pr(mo)e-fe^ (l + (/ sky Ag p 6 (m )) (24) 

is the evidence which ensures Pr(mo|d, / sky ) is normalized to unity. Recall that N s is included in the vector of 
parameters mo. 

We now move on to discuss the general case, where there are N b blobs, labelled b\ . . . bpf e , identified as containing a 
candidate source. We assume that the blobs in question do not overlap. For N T = Nb + 1 regions on the sky, making 
the approximation that the template likelihood is small when evaluated outside a blob and neglecting correlations 
between blobs, we obtain for the general case 

f 0, if N s > N h , 

Pr(d|Ag = \ (25) 
{ P r (d\N s = 0)E5> 2 ,...,6 Ws =i A blb >~ bN - nf=iP6 s (m ), if N s < N h . 

The quantity A &lb2 - bjv = is one when all indices take distinct values and zero otherwise: the sum hence generates 
all permutations of N s sources located in Nt, blobs, assuming no more than one source per blob. If there are fewer 
blobs on the sky than proposed sources, then the likelihood is very small: this would involve fitting more than one 
template within a single blob, and incurring the penalisation previously discussed. If there are at least as many blobs 
as proposed sources, then the likelihood takes the form of a sum that includes every possible association of the N s 
sources with the Nb blobs, provided that no two sources are matched to the same blob. 

Substituting Eq. [25] into Eq. |6j the expression for the approximation to the full posterior is given by 

Pr(m |d,/ Sky ) * Pr(mo)Pr(d|iV s = 0) e . /5kvJv5 g g A M,..*. J^ (mo) . (26) 

Eq. [26] is the main result of this calculation, from which all following results are derived. In the limit of a single 



isolated observation, Eq. 26 reproduces the Bayesian source detection formalism developed in Refs. [TB1IT§] . 



III. SOURCES 



In this paper, we consider two theories that give rise to localized sources in the CMB: cosmic bubble collisions 
in the eternal inflation scenario and cosmic textures. For bubble collisions, the only global parameter is N s , and 
the final result of the analysis is a one-dimensional posterior probability distribution. The first analysis of cosmic 
bubble collisions using a variant of the approximation scheme outlined in the previous section was presented in 
Refs. [HUTU], where, in addition to the location, three local model parameters (size, edge discontinuity, and amplitude) 
were included. Cosmic textures have two global parameters: N s and a measure of the symmetry breaking scale e, 
and therefore the final product is a two-dimensional posterior probability distribution. An analysis of textures, which 
also used a variant of the approximation scheme outlined above, was presented in Ref. In this study a model 

of textures with one local parameter (size) in addition to the position was considered. Previous work on testing for 
the signature of textures in the CMB was presented in Refs. [2"0"l - f2~3] . Below, we outline our models of these two 
types of sources, including the prior probability distribution on the local model parameters. For bubble collisions, we 
update the model assumptions of Refs. [9] [10] in light of improved theoretical understanding, and remove the edge 
discontinuity parameter from our analysis. 



A. Bubble collisions 



For an overview of the theory of eternal inflation and the observable effects of bubble collisions, we refer the reader 
to the reviews [BJUS]. For a detailed discussion of the expected signature of bubble collisions in the CMB, we refer 
the reader to Refs. [9, 2"6"M29| : here, we provide only a brief overview. 
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Based on the symmetry of the bubble collision spacetime, the existence of a causal boundary splitting the bubble 
interior into regions affected and not affected by a collision event, and the fact that a bubble collision is a pre- 
inflationary relic, the most general template for the temperature fluctuation caused by a single bubble collision is 
given by IIS] 

4-1 n J.\ ( z crit — ^0 COS cr it Z — Z crit „\ „,„ „, , > 

t(0, 0) = — - a + - a — cos (9 6((9 crlt - 9), (27) 

V 1 - cos 9 ciit 1 - cos crit J 

where CT it is the angular scale of the source, corresponding to the causal boundary of the collision event, zq and z cr it 
are the amplitudes at the center and edge of the template, O is the Heaviside step function, and we have centered 
the template on the Galactic north pole. In the limit of small amplitude, this is an additive contribution to the CMB 
temperature anisotropies as in Eq. [T] Theoretical work [2"B1 |2"9"] which appeared subsequent to the previous analysis 
suggests that there is no discontinuity in temperature at the causal boundary, and we therefore restrict our attention 
to z crit = 0. 

The bubble collision model contains only one global parameter, the expected number of detectable sources N s . 
This is partially a function of the properties of the potential sourcing inflation, and as such is impossible to predict 
without a model for the potential. In the context of an inflationary landscape, N s can be considered as a continuous 
parameter with some prior distribution reflecting the typical vacua produced, but without a measure for the landscape 
it is difficult to estimate even an order of magnitude for this number. It may even be quite likely that N s ^> 1, in 
which case the approximation of looking for widely separated sources is not valid; see Ref. [3D] for a discussion of the 
observational signatures in this case. In the absence of a detailed theory, we assume that N s can take values of 0(1), 
and set N s to be uniform in the range < N s < 10; as we shall see, this parameter is constrained by data, and the 
precise choice of upper limit has no effect on the analysis. In this model, ACDM corresponds to N s = 0. 

The local parameters are the collision signature's central amplitude, zo, size, C rit, and location, {#o,0o}. The 
modulations are equally likely to be hot or cold and are isotropically distributed across the sky. Theory does not fix 
the expected amplitude of the collisions, so we assume that the amplitude is uniform in the range — 1CP 4 < zq < 10~ 4 
(as stronger collisions would have been obvious in previous CMB data). Neglecting the back-reaction of the collision 
on the geometry of the bubble interior, the distribution of source sizes is proportional to sin# cr it [SI 1251 [S3 122] ■ Further 
assuming no correlation between the various local parameters, the final normalized prior on the local parameters is 

Pr(rm) = Pr(^ ) Pr(# , <f> ) Pr(0 crit ) = l — ( ( ^ ( 28 ) 

where < 9q < ir, < (j>o < 2tt, and the extrema of the size distribution, {6*™J", 6*™f t x } are chosen such that the 
bubble collisions are detectable. The lower limit, 6*"^" = 2°, stems from the fact that the CMB contains considerable 
power on the degree-scale, greatly increasing the difficulty of detection. The observable signature of a bubble collision 
with # cr it larger than the upper limit, #™^ x = 90°, would be indistinguishable from the signature of a collision of size 
180° — C rit- 



B. Cosmic textures 



For a detailed discussion of the production, evolution, and observational signature of cosmic textures we refer the 
reader to the original literature [TJ-fSl H3] . In brief, CMB photons passing through a collapsing or exploding texture 
will be red- or blue-shifted, producing an azimuthally symmetric feature on the CMB sky of angular size # cr it with 
temperature profile of the form 

t(t),4,) = -J ^ } "' »<».; t(0, « = i= ^f 1 exp f - { "lZ S ' ) ), «>«,; (29) 

where the amplitude e = 8ir 2 Gr/ 2 depends on the scale of symmetry breaking, 77, 9 cri t is the size of the texture, 
0* = \/36* cr it/2, p = {0, 1} and the template is centered at the Galactic north pole. 

Assuming that all textures are produced in a single symmetry-breaking phase transition, the texture model has two 
global parameters: the dimensionless symmetry-breaking scale e and the expected number of detectable sources A^ s . 
Texture unwinding events produce features that all have the same amplitude on the CMB sky. The total expected 
number of unwinding events depends on the particular model giving rise to textures. Simulations 23 of SU(2) 
textures indicate that we can expect to have causal access to ~ 7 textures with # cr j t > 2° in the CMB; the precise 
number of detectable unwinding events further depends on our particular realization of the background CMB, the 
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dominant source of noise in the analysis. We adopt a uniform prior on N s between < N s < 10 to encode our 
ignorance of the precise theory and the effect of our CMB realization on detectability. Again, we shall see that this 
parameter is constrained by data, and thus the choice of upper limit has no effect on the analysis. Requiring that the 
symmetry-breaking scale for textures is below the scale of cosmological inflation (bounded to be lower than ~ 10 16 
GeV by the absence of observed B-mode polarization), we place an upper bound on e of I0~ 4 . We assume a uniform 
prior on e down to 2.5 x f 0~ 5 , which corresponds to the estimated detection limit of our pipeline (details of which can 



be found in Sec. IV). Under the assumption of a uniform prior, the posterior is simply proportional to the likelihood; 
results for a different prior can be obtained easily by re-weighting the posterior. 

The local parameters for textures are the size, 6> cr i t , location, {9o, 4>o} and p, which specifies whether the texture is 
hot or cold. Textures are expected to be isotropically distributed over the sky, with a distribution of sizes proportional 
to 0~j t (see Ref. [53] for a derivation), and so the normalized prior on the local texture parameters is 



Pr(mi) = ^( _J__ 1 ] , (30) 
M^tV(2° a (50°) 2 ' 



where < 8 < n, < <fio < 2ir, and we take 2° < (9 cr ; t < 50°. The lower limit on 8 CTit results from the large power 
on degree scales in the CMB; the upper limit stems from the fact that templates with 6* cr it > 50° are large enough to 



cover the whole sky and overlap themselves, rendering Eq. 29 invalid. To marginalize over the different signs, we can 



simply sum the likelihoods evaluated at the same magnitude of e, but where the template takes opposite sign. 

We now describe our implementation of the formalism derived in Sec. [XXJ. There are two main steps in the algorithm 
- candidate-source detection and patch-based evidence calculation - in which the non-zero regions of the likelihood 
function are first estimated and then evaluated. In the first step, optimal filters matched to the signals of interest 
are used to identity the positions and sizes of the most likely sources in the dataset. We describe this procedure in 
Sec.|IV|below. In the second step, the patch-based evidence ratios, Eq.[f8j are calculated for each candidate using the 



nested sampler MultiNest [33l|3l], before combining to form the posterior, Eq. 26 on the global model parameters. 
We describe this part of the calculation in Sec. [VI 



IV. CANDIDATE DETECTION WITH OPTIMAL FILTERS 



In order to effectively approximate the full posterior distribution describing the population of candidate sources, it 
is necessary to first locate the regions that provide significant contributions to the likelihood. We follow the approach 
of Ref. [H] , using optimal filters defined on the sphere that are matched to the source profile of either bubble collision 
or texture signatures. First, we construct optimal matched filters for the purpose of detecting candidate sources 
embedded in full-sky WMAP 7-year data and assess their performance. Second, we briefly describe the optimal- filter- 
based algorithm for detecting sources of unknown and differing sizes, highlighting differences between the bubble 
collision and texture cases. Third, we calibrate the algorithm on an end-to-end simulation of WMAP observations, 
before assessing its sensitivity. 



A. Optimal bubble collision and cosmic texture filters 



We construct two sets of matched filters: one set that enhances the contributions of bubble collision signatures and 
one set that enhances the contributions of texture signatures. The matched filters are constructed to enhance the source 
profile in a specified stochastic background. A stochastic background of CMB fluctuations is assumed, characterized by 
the lensed ACDM power spectrum that best fits the WMAP 7-year data, baryon acoustic oscillations and supernovae 
observations [35] (hereafter referred to as the lensed WMAP7+BAO+H0 power spectrum). The bubble collision and 
texture source profiles for which we search are relatively large-scale; thus we consider spherical harmonics up to the 
band-limit i max — 256 only. Since we eventually apply these filters to W-band WMAP observations, we assume 
observations are made in the presence of a Gaussian beam of full-width-half-maximum FWHM = 13.2 arcminutes 
and isotropic white noise of Nt — 0.02 //K 2 . The optimal matched filters are constructed in harmonic space: thus the 
assumption of an azimuthally-symmetric beam profile and isotropic noise simplifies their construction and application 
considerably, while remaining highly accurate for the relatively low band-limit considered. Once the source profile 
and stochastic background are defined, the filters are constructed on the sphere as outlined by Ref. [55] , 

In Fig. [2] and Fig. [3] we show the matched filters recovered for the bubble collision and texture profiles, respectively, 
for a range of source sizes. Notice that the bubble collision filters on smaller scales contain a central broad hot region 
to enhance the main bubble collision contribution, surrounded by hot and cold rings to enhance the transition from 
the collision to the background. On larger scales, however, the matched filters contain only the hot and cold rings 
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(a)0 crit = 5° (b)0 crit = 20° (c)0 crit = 60° 

FIG. 2. Matched filters optimized to bubble collision signatures of varying size embedded in a ACDM CMB background. 




(a)0 crit = 5° (b)flcrit = 10° (c)e crit = 20° 



FIG. 3. Matched filters optimized to cosmic texture signatures of varying size embedded in a ACDM CMB background. 

that enhance the transition. Since the CMB has more power on large scales, the matched filters on large scales do not 
respond to the large-scale features of the bubble collision signature but rather the transition region near the location 
where the template goes to zero. The texture source profile has a smooth, Gaussian transition to the background, 
and consequently the matched filters recovered for textures contain only a central broad region, without any strong 
contribution from the perimeter of the profile. 

The matched filters constructed are optimal in the sense that no other filter can yield a greater enhancement of 
the signal-to-noise ratio (SNR) of the filtered field. It is possible to calculate analytically the SNR of the filtered field 
for various filter types, as derived for example by Ref. [12] . In Fig. [4] we plot the SNR computed for bubble collision 
and texture profiles for the unfiltered field, for the optimal filters constructed here, and for needlets [37l|38], which 
have been used previously to detect candidate sources [SHTT]. Note that the lack of a sharp transition in the texture 
template means that the SNR for textures are lower than those of bubble collisions. Nevertheless, it is clear that 
matched filters yield the highest SNR, in accordance with expectations. 



B. Candidate object detection algorithm 



Although we have constructed optimal filters for a range of source sizes, we have not yet addressed the problem 
of detecting sources of unknown and differing sizes. We adopt the algorithm described in detail by Ref. [12] for this 

purpose, which we review here briefly. First, matched filters are constructed for a grid of source sizes R € {0crit}fe=T" ■ 
All filters are then applied to the full-sky observed data by convolving the matched filter kernel with the observed 
data, which may be computed efficiently in harmonic space (see e.g. Ref. [H]). Significance maps are then computed 
by normalizing the filtered field to the mean and standard deviation of filtered fields computed from realizations of 
the background process (i.e. CMB fluctuations and instrumental noise) in the absence of sources. The significance 
maps are then thresholded (the calibration of threshold levels is discussed below), before potential candidate sources 
are found from the localized peaks of the thresholded significance maps. Potential candidate sources are eliminated 
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if a stronger source is found on adjacent scales, where the set of scales adjacent to scale R is defined by the set 

{-Radj € {#crit}fc=i rit : l-^adj — R\ < ^adj}, i-e. where the distance between R and i? a dj is less than the parameter 6* a dj. 
Once candidate sources are detected, the parameters of the source size, location and amplitude are estimated from 
the corresponding filter scale, peak position of the thresholded significance map and amplitude of the filtered field, 
respectively. 



The construction of optimal filters is implemented in the S2FIL code [36] (which in turn relies on the codes S2 [39 

and HEALPix [10]), while the COMB code has been used to simulate bubble collision signatures embedded in a CMB 
background^] The candidate object detection algorithm described here is implemented in a modified version of S2FIL 
that will soon be made publicly available. 

There is no guarantee that the peak in the filtered field across scales will coincide with the size of the unknown 
source. Nevertheless, for bubble collision signatures embedded in the CMB this has indeed been found to be the case 
[12] . For texture profiles, however, we have found this phenomenon to hold below scales 8 clit ~ 10° only. Through 
numerical experiments we found that the difference in the behavior of the filtered field between bubble collisions and 
textures on large scales is due to the absence of a well-defined transition region from the source to the background in 
the texture profile. For large texture sizes, not only is there no peak in the filtered field at the scale of the unknown 
source, but the SNR of the filtered field does not drop off rapidly when applied to nearby scales (Fig. El), as is the 
case for bubble collisions. We trivially modify the candidate detection algorithm described above to account for this 
behavior. For textures sizes below cr it = 10° we look across adjacent scales as usual to find the most significant 
potential candidate source, whereas for sizes above 9 CIlt = 10° we do not (by a judicious choice of the adjacency 
parameter 8^3 there is in fact no need to modify the algorithm, as described below). 

Although the candidate detection algorithm considers a grid of candidate scales R, it is overwhelmingly probable 
that the signal for any given source peaks at scales between the samples of the grid. It is thus important to examine 
how sensitive the matched filters are to small errors in the source size. In Fig. H (b) and (d) we plot SNR curves for 
matched filters constructed on the grid of candidate scales. For bubble collision profiles, a sharp degradation in the 
SNR away from the scale used to construct each filter is clearly apparent. For texture profiles the degradation is much 
less pronounced, especially at large scales (as discussed above). Provided that the # cr it grid is sampled sufficiently 
densely, the matched filters for both bubble collision and texture profiles remain effective and are superior to needlets. 



We define the parameters of the optimal-filter-based candidate object detection algorithm here and calibrate the 
threshold levels for WMAP observations. Throughout the calibration we apply the WMAP KQ85yr7mask [3T], since 
this is the mask adopted when analysing WMAP data subsequently. We select the less-conservative KQ85 mask so as 
to reduce the variance in reconstructing large-scale information masked by the Galactic sky cut [32]. This choice has 
the additional advantage of revealing more of the sky and hence increasing the number of candidates. The Baycsian 
stage of the pipeline will assess the overall evidence for the source model in each candidate region: if any candidate is 
found to contribute evidence in favor of the source model, then it will be examined closely for frequency-dependence 
signifying potential foreground contamination. 

For bubble collisions we consider the grid of scales set out in Ref. [12] . as defined in Table [I] (left). The 6* cr ; t prior 
range is smaller for textures than for bubble collisions - the texture profile extends well beyond # cr it , covering the full 
sky for 8 alt > 50° - and hence for textures we consider a smaller grid of scales, also defined in Table [I] (right). Since 
we found the matched filters for textures to be sensitive to a large range of nearby scales, we nevertheless remain 
sensitive to the full prior range of sizes. The SNR curves for the matched filters constructed for these scales for bubble 
collisions and textures are shown in Fig. [4] (b) and (d). These grids of scales are thus sufficiently sampled to ensure 
that the matched filters remain effective for scales between the samples of the grid. We set the adjacency parameter 
to # a dj = 5° for both bubble collisions and textures. For textures this ensures that we look across scales for sizes 
below 9 cl a = 10° but not above, whereas for bubble collisions we always look across adjacent scales. 

We use 3,000 Gaussian CMB simulations to calculate the mean and standard deviation of the filtered field at each 
scale in the absence of sources, in order to compute significance maps. For these simulations, and for the WMAP 
data analysed subsequently, we perform Wiener filtering to recover spherical harmonic coefficients with I < 10 from 
masked CMB maps [32], where we adopt a Gaussian prior for the harmonic coefficients specified by the lensed 
WMAP7+BAO+H0 power spectrum. Note that this differs from the maximum likelihood reconstruction [43] of 
harmonic coefficients performed by Refs. [5HTT] and the cut-sky estimation performed by Ref. [T2] . This can alter the 





C. Candidate object detection calibrated to WMAP 



3 S2FIL, S2 and COMB are available from 



http : //www . jasonmcewen . org/ 



while HEALPix is available from http://healpix.jpl.nasa.gov/ 
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FIG. 4. SNRs of bubble collision (top row) and texture (bottom row) signatures of varying size with amplitude 100 fiK 
embedded in a ACDM CMB background. SNR curves are plotted for matched filters (solid blue curve), needlets [37] with 
scaling parameter B = 1.8 for a range of needlet scales j (dot-dashed black curves) and for the unfiltered field (dashed red 
curve). In panels (b) and (d) SNR curves for the matched filters constructed at a given scale and applied at all other scales 
are also shown (thin solid blue curves). The scale for which the filters are constructed may be read off the plot from the 
intersection of the heavy and light solid blue curves. Provided the # cr it grid is sampled sufficiently densely, the matched filters 
remain superior to needlets. 



spherical harmonic coefficients recovered on large scales, and thus the detected candidate sources, in a non-negligible 
manner. However, Wiener filtering should give the most reliable reconstruction of the large-scale harmonic modes 
02]. 

Finally, we calibrate the threshold levels N aR applied to the significance maps for each filter scale from a realistic 
WMAP simulation that does not contain embedded sources. The thresholds are chosen to allow a manageable 
number of false detections while remaining sensitive to weak sources. For this calibration we use a complete end-to- 
end simulation of the WMAP experiment provided by the WMAP Science Team [31] . The temperature maps in this 
simulation are produced from a simulated time-ordered data stream, which is processed using the same algorithm as 
the actual data. The data for each frequency band are obtained separately using simulated diffuse Galactic foregrounds 
and CMB fluctuations, and include realistic noise, smearing from finite integration time, finite beam size, and other 
instrumental effects. We use the foreground-reduced W-band simulation for calibration. The threshold levels N CTR 
are selected to allow at most three false detections on each scale on this simulated map (recall that detections on one 
scale can be eliminated by stronger detections made on adjacent scales). These threshold levels are less conservative 
that those set by Refs. [5rTT2"]. in order to increase the number of false detections passed by the candidate source 
detection stage of the analysis pipeline and hence improve its sensitivity. The calibrated threshold levels for both 
bubble collisions and textures are shown in Table [I] Once candidate sources are detected by the optimal-filter-based 
detection algorithm, we discard those objects that are significantly masked. For the WMAP end-to-end simulation, 
we make 12 false detections of candidate bubble collisions and 4 false detections of candidate textures. 
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TABLE I. The # cr it grid and threshold levels N aR adopted for the optimal-filter-based candidate source detection algorithm 
for bubble collisions (left) and textures (right). Threshold levels are calibrated to the WMAP end-to-end simulation to allow 
at most three false detections on each scale. 



In the analysis of the WMAP end-to-end simulation (and in the analysis of the WMAP data considered subse- 
quently) , some of the bubble collision candidates that we detect differ from those found with optimal filters previously 
[12] . This is expected since we now use Wiener filtering to recover spherical harmonic coefficients, have included a 
more accurate model of the WMAP noise in the optimal filter construction (noise was neglected previously), and have 
reduced the threshold levels in order to increase the sensitivity of the entire pipeline. Further, the thresholding-based 
nature of the candidate source detection algorithm means that small differences in the filtered field can have an 
impact of the final candidates detected if regions move below or above the threshold. In the cases where candidates 
disappeared, we nevertheless found peaks in the filtered field; these were simply no longer above the threshold or 
were eliminated by stronger detections on nearby scales or positions. Given these differences are due to improvements 
made to the pipeline, the results given here are to be preferred to those presented previously. 

Following the candidate source-detection stage of the analysis pipeline, we pass to the Bayesian stage an estimate of 
the domain of parameter space over which the likelihood is expected to be non-negligible. These regions of parameter 
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FIG. 5. Sensitivity of the optimal-filter-based candidate detection algorithm, with completeness curves plotted for one, two and 
three standard deviations from the 50% completeness curve. The completeness curves are computed in the following manner. 
For each scale # cr it we compute the source amplitude (zo for bubbles; e for textures) that would be required to ensure that 
the SNR reaches the threshold specified in Table [I] This level defines the 50% completeness curve, since we expect half of the 
sources with this amplitude to fall below the curve and half to fall above. Similarly, we compute approximate completeness 
curves for one-, two- and three-standard-deviation differences from the 50% completeness curve. The probabilities quoted on 
each completeness curve are computed assuming a Gaussian distribution of the filtered field at the source position. 



space are estimated from each of the candidate sources detected. For the size of each source, the relevant region is 
determined first by finding the range of nearby filter scales for which significance maps exceed their threshold level 
at the source position. This range of scales is extended to the next smallest and largest filter scales (or the edge of 
the prior if encountered) to yield an estimate of the range of scales over which the likelihood is non-negligible. For 
example, for a bubble collision candidate found to be significant by the 40° and 45° optimal filters, we would estimate 
the likelihood to be non-negligible over the range 35 - 50°. To estimate the integration limits of the central positions 
for each source, we first find all pixels within a radius r of the source position estimated by the optimal filters, where 
r is 25% of the maximum source size estimated from the previous step. The extrema {#™ m , 6*™ ax , </> nm , 0™ ax } of 
these pixels are found, and the source positions are then sampled from the region defined by #™ ln < 9q < 6* nax , 
0O 1111 < 0o < <^o Qax - Tests of these assumptions are included in the suite of pipeline tests detailed in later sections. 

We conclude this section by assessing the level to which the optimal-filter-based candidate detection algorithm is 
sensitive for each source type. In previous studies, simulations were performed for this purpose [9l I10[ I12j . Here 
we instead take a probabilistic approach based on the analytic SNRs of the filters computed previously (see Fig. [4j. 
This allows us to probe the source size-amplitude parameter space at higher resolution and accuracy than would be 
achievable with modest numbers of simulations (to reach an equivalent resolution and accuracy through simulations 
would be extremely computationally demanding). In Fig. [5] we plot the sensitivity of the matched filters constructed 
for bubble collisions and cosmic textures. These plots are produced as follows. For each scale cri t we compute the 
source amplitude (z for bubbles; e for textures) that would be required to ensure that the SNR reaches the threshold 
specified in Table |TJ This level defines the 50% completeness curves shown in Fig. [5] since, in the presence of noise, we 
expect half of the sources with this amplitude to be detected and half to be missed. Similarly, we compute approximate 
completeness curves for one-, two- and three-standard-deviation differences from the 50% completeness curve (note 
that the probabilities quoted on each curve are computed assuming a Gaussian distribution of the filtered field at the 
source position). For the 50% completeness curve, the bubble collision matched filters are sensitive to zq ~ 10 -44 , 
while the cosmic texture matched filters are sensitive to £ ~ 10~ 42 . Note that the sensitivities computed in this 
manner are similar to those computed previously through simulations [HI [101 both in terms of the sensitivity 
levels obtained and the shape of the sensitivity regions. Further, we see that optimal matched filters are ~ 1.7 times 
more sensitive than needlets for detecting bubble collision signatures, as found previously [9) fTTJl [12] . 



V. ADAPTIVE-RESOLUTION EVIDENCE CALCULATION 



Modern CMB experiments map the sky with extremely high resolution: the beam of the Planck experiment in 
the main CMB bands is expected to be ~ 5' [17], resulting in maps pixelated on the arcminute scale. While this 
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1 10 100 

patch radius (degrees) 

FIG. 6. The memory needed to calculate the covariance matrix for patches of a given size at HEALPix resolutions ranging from 
A^sidc = 8 (lightest-blue, lowest curve) to -/V s id e = 2048 (highest, darkest curve). The effects of pixel size are visible at small 
patch radii and low resolutions. 

is necessary for pinning down the secondary CMB anisotropics at small scales, it means that calculating pixel-space 
covariance matrices becomes extremely memory-intensive. We illustrate this point in Fig. [6j which shows the memory 
needed to calculate covariance matrices from 1° to 180° in radius at HEALPix resolutions ranging from iV S id e = 8 
to Nsido = 2048 (i.e., Planck resolution) Q It is clear that the memory costs, which to a good approximation rise as 
angular radius to the fourth power, make processing even relatively small patches prohibitive at full Planck resolution. 

In previous work \Q UOj , we chose to truncate both our patches and the integration limits of # cr it to the maximum 
radius invertible with our memory constraints. While this allowed us to at least partially process almost all features at 
full WMAP resolution, it meant that we were unable to probe the large-0 cr ; t region of parameter space for which the 
prior for bubble collisions is highest. If we are to do so, it is clear that the input maps must be processed at degraded 
resolution: the larger the patch, the lower the resolution at which it can be processed. The maximum memory 
accessible per core in this analysis is ~ 90 GB, which means that the full-sky covariance matrix can be inverted at 
iVgide = 64: this is therefore the minimum resolution at which any feature will be processed. The degradation process 
will now be described in detail, along with the suite of tests performed to assess its performance. 

A. Processing maps 

The candidate detection stage returns estimates for the size and position of features of interest within a map, defining 
the set of patches to be considered in the evidence calculation. The memory cost of computing each covariance matrix 
is derived from the number of unmasked pixels within the patch: if this is greater than the memory available, the 
patch must be processed at a degraded resolution. 

In the HEALPix pixelization scheme, each step down in resolution reduces the pixel count, N p i x , by a factor of four. 



4 The quantity plotted corresponds to a total number of matrix elements equal to ~ 1.5 X N^ ix - Our algorithm calculates the full N^ ix 
covariance matrix in order to make use of the LAPACK inversion routines |44| . then compresses the inverted matrix to a 1-D array 
containing its upper triangle to reduce memory costs while sampling. 
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■A^i side,deg 




Pixel Scale (arcmin) 


Smoothing FWHM 


256 


786432 


13.7 


33.9 


128 


196608 


27.5 


67.7 


64 


49152 


55.0 


135.2 



TABLE II. The full-widths at half-maximum of the Gaussian kernels used to smooth input maps prior to degradation. Also 
tabulated are the pixel count at each resolution, and an approximate pixel scale, derived assuming each pixel is square. 



As the covariance matrices are A^ ix in size, an estimate of the number of steps down required is therefore given by 

/log(m ost /m max )\ 
n deg = ceiling — , (31) 



log 16 

where m ost and m max are the estimate of the memory required and the memory available. The corresponding estimate 
of the finest resolution at which the patch can be processed is then 

. r -^side.full / Qr) x 

A/ Si dc,dcg - 2 „ dcg ■ [M) 

This estimate is tested by re-counting the number of unmasked pixels in the patch at the new resolution. An 
-^side,deg mask is created by averaging within degraded-resolution pixels: any N s - lt i e ^ eg pixel which is more than half 
masked at the input resolution is considered to be masked. A precise calculation of the memory cost of the degraded 
covariance matrix is then made: if this is below the memory threshold, as expected, the degraded resolution is accepted 
and the algorithm proceeds; if not, the resolution is decreased once more. 

Once the required resolution has been determined, the input CMB temperature map can be degraded. As the 
CMB is a smooth field, it is not sufficient to simply average within Af s i<j e ,deg pixels: doing so will introduce large 
pixelization effects, which will act as an extra noise term unaccounted for in the pixel-pixel covariance matrix. This 
can be avoided by smoothing such that the input map is smooth on the degraded-pixel scale prior to reducing the 
resolution. This is equivalent to introducing a band-limit, f max , in harmonic space. Choosing the band-limit is a 
balance. If the smoothing scale is set too large (i.e., £ max is too low), too much information will be discarded with 
each degrade and performance will suffer. If the smoothing scale is set too small (i.e., smaller than an A^ide.deg pixel), 
the smoothed maps will contain pixelization artefacts. 

The choice of £ max is somewhat arbitrary, but experimentation shows that the degradation is stable if the harmonic- 
space Gaussian smoothing kernel is 1% of its maximum at £ ma ,x = 2-/V S idc,dcg- This defines a smoothing scale at each 
resolution: the FWHM, / d eg> of the pixel-space kernel is given by 

_ / 8 log 2 log 100 

; dcg ~ V I U +TT ' { ' 

y ^max ^max i ^ ) 

Assuming the 12AT s 2 ide pixels in a HEALPix map of given resolution are flat and square (a safe assumption at high- 
resolution), the smoothing scale is approximately 2.5 times the size of a pixel, and the maps are clearly smooth on the 
pixel-scale. The smoothing kernel sizes used for the three degraded resolutions considered in this work are listed in 
Table [Tlj Note that for speed and simplicity the smoothing is carried out on the full-sky, rather than within a patch: 
the time taken to smooth in pixel-space scales poorly with patch size. 



B. Calculating the degraded evidence 



Care must be taken when calculating the covariance matrix for use in the likelihood: the covariance matrix must 
include as faithful as possible a representation of the components of the data. We must therefore capture every 
important feature of both the CMB "signal" (the CMB is, in fact, the dominant noise in the analysis) and instrumental 
noise measured in the input map, as well as the effects of the degradation process. It is helpful to break up the full 
covariance matrix, C, into its CMB, S, and noise, N, constituents, as they are morphologically different. 
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FIG. 7. Window functions of the various smoothing kernels appearing in the adaptive-resolution analysis. Shown are the 
WMAP7 W-band beam (approximated as a Gaussian, solid line), the Gaussian smoothing applied before degradation to 
A^side.dog = 128 (dashed), and the pixel window function at this resolution (dotted). Note that the pixel window function is 
only defined up to I = 4iV s id ej d eg , the maximum multipole allowed by the HEALPix software. 



1. CMB covariance 



The CMB, as a correlated random field on the sphere, is most simply defined in harmonic space by its power 
spectrum, Cg. At full resolution, the CMB power is smoothed by the instrumental beam, which we approximate 
in this analysis with a Gaussian of FWHM /wmap- In patches that are processed at reduced resolution, the CMB 
signal is also smoothed by the anti-aliasing beam (another Gaussian, of FWHM /deg) and further by the pixel window 
function of the degraded resolution map. This final effect, shown for the relevant beams in Fig. [7J is small but 
non-negligible. Ignoring the degraded pixel window function means that the covariance contains an overestimate of 
the CMB power - in our analysis, a noise term - and log-evidences can be underestimated by as much as 1 when 
degrading to iV si d e ,deg = 64. 

Taking all of these effects into account, the CMB covariance between two pixels i and j is 

S v = E 2 ^ C ^( cos d »)Bj, wmap (34) 

t 

if the patch is processed at full resolution, and 

Sij = E -^^ CeP ?( COs9i j) B twMAP B tdc e W £,dcg ( 35 ) 

e 

if the patch must be processed at degraded resolution. Here, Pt(x) are the Legendre polynomials, W^,deg is the iV s ide,deg 
pixel window function, and the Bg are the WMAP and anti-aliasing Gaussian beams, represented in harmonic space 
as 



16 log 2 



where / is the relevant FWHM. 



2. Instrumental noise covariance 



At full resolution, the WMAP noise is uncorrelated in pixel space (see, e.g., Ref. [IS]). The noise covariance matrix 
is therefore diagonal, and can be written as 



(37) 
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where ctwmap = 6.549 mK is the RMS noise of the W-band detectors, iV bs.j is the number of times pixel i has been 
observed, and 5y is the Kronecker delta. 

The anti-aliasing smoothing applied as part of the degradation process induces correlations between the noise 
measured in each pixel. Coupled with the variations in sky coverage represented by the N a i, s values, this makes the 
exact pixel-space covariance more difficult to write down. Progress can be made by separating the full-resolution noise 
model into isotropic white noise, calculated from the power spectrum in harmonic space, which is modulated by a 
map encoding the variations in sky coverage, which belongs entirely in pixel space. The effects of the smoothing and 
degradation can then be applied individually to each component, and recombined in the covariance matrix. 

Uncorrelated pixel noise is identical to the noise generated by white (i.e., flat) noise in harmonic space. To produce 
the correct RMS noise, the amplitude of the white noise power spectrum, Ne, should be set such that 

N e = g WMAP^ P iX ; (38) 
Nobs 

where f2 p i X is the pixel area, and we have incorporated the mean number of observations per pixel, iV b s , into the 
definition. After smoothing and degradation, the noise power spectrum is no longer flat, having been multiplied by 
both the anti-aliasing beam and the degraded pixel window function (both squared). 

Having absorbed the mean number of pixel observations into the isotropic noise, to include the effects of varying 
sky coverage we need only consider relative changes in the number of observations in each pixel. These are captured 
by generating a map of \J iV bs/^obs,i- After degradation, this map will have been smoothed and degraded exactly as 
the input temperature map; its mean at any resolution is one. The noise covariance matrix at degraded resolution is 
then 

Nij = /A Nohs M £ 2 ^N e P i (co S e^Bl^Wl^, (39) 

V iV obs,i-'V b S j e 47r 

where the separation into components residing in pixel and harmonic space is clear. Note that this expression does 
not include the instrumental beam, as this is detector noise. 



3. Full covariance 



The sum over £ in the covariance matrix calculations is strictly a sum to infinity, but is approximated with a sum 
from < £ < ^max.cov, truncated at a point where all significant contributions have been included. At full resolution, 
the data are noise-dominated at high-£ (see Fig.ls]), so the covariance is not particularly sensitive to the precise value 
of 4nax,cov) provided it is larger than ~ 1000. We use f max ,cov = 10240 After degradation, the CMB and noise power 
spectra have been additionally damped by the extra smoothing and pixel window function. These effects combine to 
reduce the map power to ~ 0.7% of its maximum at £ max = 2iV s idc.dcg7 but the kernel's long tail means that there 
is still information at higher-£. If this is ignored, the maps contain more information than is accounted for in the 
covariance matrix. The signal-to-noise ratio of any relic present is therefore artificially boosted, and the evidence is 
over-estimated: truncating at 2./V s jde,deg yields over-estimates of ~ 5 in the log-evidence at iV s ide,deg = 64. Tests reveal 
that evidence ratios are stable provided multipoles £ > 3N s id e ,deg are included in the covariance matrix. We choose 
to be conservative and push this to £ m a,x,cov = 4A s id c ,dcg- 

The two components of the covariance are added to form the full covariance matrix. 

C = S + N. (40) 

This matrix must be inverted for use in the likelihood function; the inversion is carried out using the Cholesky 
decomposition implemented in the LAPACK software package. The full-resolution covariance matrix has a prominent 
diagonal due to the uncorrelated pixel noise, and hence is readily invertable. However, at lower resolutions the 
additional smoothing can make the covariance matrix nearly singular. Diagonal regularizing noise is therefore added 
at the level of 2 fiK 2 to all degraded-resolution covariance matrices to aid inversion: we have checked that this does 
not affect the results of the likelihood calculation. 



5 The WMAP7 observations |35l extend the power spectrum measurements to I = 1200, but the CMB signal-to-noise ratio for the 
W-band is below one for I > 600. Setting £ maXjCO¥ = 1024 ensures the CMB contribution is characterized well into the noise-dominated 
regime. 
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FIG. 8. Power spectra of the WMAP7 best-fit CMB signal (solid lines), WMAP7 noise (dotted) and a single simulated bubble 
collision (dashed) with amplitude zq — 5 x 10 -5 and angular radius # cr it = 7°. The power spectra are plotted at full WMAP 
resolution, iV s idc = 512 (left), and after smoothing and degradation to iV s id e ,deg = 128 (right). 



C. Calculation of patch-based evidence ratios with MultiNest 



To form the full posterior (Eq. |26[), we must calculate the patch-based evidence ratios (Eq. 18) for each blob. Our 
pipeline uses the MultiNest sampler [331 [M], which outputs the normalized posterior and the evidence for the data 
in each patch under a specific model. This is not precisely the information required in Eq. |18[ which integrates the 
product of the likelihood and the full-sky prior. However, we can use Bayes' theorem in the patch to convert the 
output posterior into a likelihood according to 

p / I j\ Pr(m )Pr fc (mi)Pr b (d|m , mi ) 

Pr 6 (m ,mi|d) = Pr b (d) ' ' ' 

Here, the subscript b indicates that the probabilities are formed using only the data for a single blob. In particular, 
we identify 

Pr b (d|m , mi ) oc g-^C^"'^-^ (42) 
as the quantity necessary for the patch-based evidence ratio. The normalization assumed in Eq. [41] is 

dm dm 1 Pr(m )Pr b (m 1 ) = 1, (43) 



which implies that Pr b (mi) is related to the prior Pr(mi) in Eq. 18 by an overall normalization: 

F b = ^ r(mi) . = 1 - /"dm dm 1 Pr(mo)Pr(m 1 ). (44) 
rr b (mij J b 



We also have 



Pr b (d) = / dm dm 1 Pr(m )Pr b (m 1 )Pr b (d|m ,m 1 ). (45) 

Jb 

Using these expressions, we can solve for the patch-based evidence in terms of known quantities: 

" 6(mo) = Fb X Pr^O) X iMmo) / dmiPrb(m °' mi|d) ' (46) 

where Pr b (d|0, 0) is the likelihood in the patch with no sources. 

All steps of the algorithm are parallelised wherever possible using OpenMP and MPI to take advantage of both 
shared- and distributed-memory clusters. 
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Bubble 6U it (°) 


-^^side,deg 






7 


512 


41618 


12.6 =b 0.3 




zoo 


1U044 


Iz.o ± U.O 


15 


256 


40715 


14.2 ±0.3 


15 


128 


10314 


14.0 ±0.3 


30 


128 


49894 


13.2 ±0.3 


30 


64 


12619 


13.1 ±0.3 



TABLE III. Tests of the stability of the degraded evidence values. Three maps, each containing a small, medium or large 
simulated bubble collision, are used to examine how the evidence ratio changes when a patch is degraded from iV s idc = 512 to 
256, 256 to 128 and 128 to 64. 



VI. ADAPTIVE-RESOLUTION TESTS 
A. Stability of degraded evidence values 

The adaptive-resolution analysis pipeline was tested thoroughly to estimate the effect the degrading has on the 
calculated evidence values. Naively, we would expect the evidence not to change (on average) for resolutions at which 
the feature is well-sampled, i.e., for which the pixel scale is much smaller than the feature itself. Once the resolution 
has decreased enough for the feature - and hence the template used to fit it - to appear pixelated, the evidence should 
begin to drop off sharply. In harmonic space this is equivalent to requiring that enough modes are left intact by the 
pre-degradation smoothing that the template power spectrum can be discerned. 

The ideal test in this situation is to create a simulation containing a sufficiently large template to be un-pixelated 
at the lowest resolution at which a feature can be processed: in this case, -/V s ido = 64. This feature should then be 
processed at all resolutions considered, from the highest (WMAP-resolution) to the lowest, to determine how the 
evidence behaves. Unfortunately, the very nature of the problem makes this is impossible: such templates would 
require enormous covariance matrices at full resolution. 

Progress can be made by breaking the test into parts. Simulations containing templates - strong bubble collision 
signatures, in this case - are generated on a range of angular scales. Each simulation is processed twice: first at 
the highest resolution possible, then again at one resolution lower. This indicates how the evidence changes when a 
patch is not processed at its "ideal" resolution, which could occur if the maximum angular size, # cr it, is overestimated. 
The lower-resolution evidence values are calculated with four times fewer pixels: if the evidences returned by the two 
runs do not differ greatly, we can be confident that small reductions in resolution do not affect the evidence values 
returned. 

Three maps are generated containing single bubble collision signatures of angular scale # cr i t equal to 7°, 15° and 
30°, each of which contains approximately the same number of pixels at = 512, 256 and 128, respectively. In 

each case the signatures have the same amplitude, z = 8 x 10~ 5 , and position, (# o ,0o) = (45°, 45°); the maps also 
contain the same CMB and noise realizations, and are plotted in Fig. [9} The evidence is calculated once at the highest 
resolution possible - A^idc of 512, 256 and 128 for the 7°, 15° and 30° collisions, respectively - and once again after 
degrading one step in resolution - i.e., for N s ^ c of 256, 128 and 64. In each pair of tests the same ranges in size and 
position are sampled, and no mask is used. The only small difference comes in the pixels included in the patch, as 
the lower-resolution patches sample a slightly larger region than the high-resolution patches]^] 



The results of the tests are shown in Table III In each case, the evidence calculated at high resolution matches the 
evidence at low resolution to within MultiNest precision. We conclude that the adaptive- resolution analysis produces 
stable evidence ratios for the resolutions considered in this work. Note that the evidence does not always decrease: in 
fact it can increase. This is because the realization noise (i.e., the combined CMB and noise signal) is different after 
degradation: the bubble collision signature is compared to larger-scale modes at lower A^idc.dog- Further, note that 
the log-evidence errors tabulated are the typical random variations due to sampling, estimated by repeatedly testing 
an individual patch with initial conditions set from different random seeds. MultiNest also provides a statistical 
estimate of the error in an evidence calculation, derived from the relative entropy of the samples (see Refs. [3"3l 146] ). 
We find these estimates (typically ~ 0.1 in log-evidence) to be subdominant to the variations due to sampling. 



When defining the patches using the HEALPix query _disc subroutine, the "inclusive" option is set to include all pixels which fall even 
partly within the radius we sample. 
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7° Bubble Collision 15° Bubble Collision 30° Bubble Collision 



1 . 



-504'— — 655 /J.K -504 — — 655 /J.K -504 — — 669 yuK 

FIG. 9. The maps used to test the adaptive-resolution analysis. Three maps are generated, each containing a single bubble 
collision of radius # cr it = 7°, 15° or 30°. The maps contain identical CMB and noise realizations, and in each case the bubble 
collision is placed at located at (6o,4>o) = (45°, 45°) with amplitude zo = 8.0 x 1CP 5 . The 7° and 30° collision maps are also 
used to test the effects of neglecting correlations with data outside the patch and the restriction of template locations to the 
regions highlighted by the candidate detection stage. All plots are shown at the same scale, and are 67° on a side. 



B. Robustness to smoothing-induced contamination 



The simulations used to test the adaptive-resolution pipeline contain only a CMB realization, a noise realization and 
a bubble collision template; no mask was used. Real CMB datasets also contain foregrounds (or foreground residuals 
after component separation), the worst of which are masked. Due to the need to smooth prior to degradation, there is 
a potential for contaminants to leak from behind the maskQ While it is possible to mitigate this effect by extending the 
mask (and potentially using a smoothing kernel that is localized in pixel space) [42] , it is highly undesirable to discard 
hard-won data. The likely scale and amplitude of any smoothing-induced contamination is therefore investigated to 
determine whether the mask should be extended. 

The masks recommended by the WMAP team |41j comprise two components: a Galactic cut (of varying conser- 
vatism) and a point-source mask. The point-source mask is created from a range of external and internal catalogs (as 
listed in Ref. |47j). and is updated with each data release. The point-sources, which appear in the data as approxi- 
mately Gaussian peaks with the same FWHM (0.22°) as the instrumental beam, are masked by excising a region of 
radius ~ 0.6° centered on each point-source. A small number of the strongest sources are more aggressively masked, 
out to a ~ 1.2° radius. The cut made to remove the extended emission of the Milky Way is much larger, and forms 
an irregular band ~ 20° — 40° in width centered on the Galactic Plane. For clarity, the effects of the two components 
of the mask are investigated individually. 

We can estimate the effects of point-source smoothing using a very simple test. A W-band point-source is simulated 
by placing a normalized 1 /iK delta-function at a position taken from the WMAP point-source catalog [H] , then 
convolving it with a Gaussian of 0.22° FWHM. Since smoothing is a linear process, this can then be scaled to 
investigate the effects of sources of different amplitudes. This map is then smoothed and degraded to each scale used 
in the analysis (i.e., iV S id e = {256,128,64}), and masked using the degraded point-source mask at each resolution. 
The resulting maps, plotted in Fig. [Toj can then be scaled to mimic a source of a given temperature. The plots show 
that, at all resolutions considered, the maximum contamination injected into a single pixel is a few thousandths of the 
point-source's amplitude. Assuming that such sources have amplitudes of 100 — 1000 /iK [35], these results suggest 
that our degradation technique induces contaminants of at most 1 — 2 fiK into fewer than 10 pixels. This level is 
completely subdominant to the CMB signal, and so should not affect the analysis. We therefore need not extend the 
point-source mask when smoothing and degrading. 

An estimate of the contaminants leaked from the Galactic cut can be obtained in a similar fashion. In place of 
the simulated point source, we must make an estimate of the Galactic foreground residuals. Modeling this precisely 
is difficult, so, following Ref. [49], 1% of the difference between the V-band signal and the WMAP7 Internal Linear 



The smoothing procedure is fastest in harmonic space, where it is a multiplication rather than a convolution. Inclusion of the mask in 
this procedure is complex, and it is simplest to perform the smoothing on full-sky data. 
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W-Band Point-Source 



1/xK Point-Source (N side = 256) 
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1/zK Point-Source (N side = 128) 1/xK Point-Source (N side = 64) 




0.0016 



FIG. 10. A typical WMAP W-band point-source contaminant (top-left, plotted using the inverted point-source mask), and the 
effects of smoothing a unit-amplitude simulated point-source and degrading to A^ide.deg = 256, 128 and 64. All low-resolution 
plots are shown with a degraded mask applied, which limits any residuals to O(fiK). All plots are 12.5° on a side. 



Combination (ILC) map [H] is used for illustrative purposes. This combination is indicative of the morphology and 
amplitude of residuals within the ILC: there are contaminants of around 50 times the amplitude visible to the eye in 
the WMAP foreground-reduced maps. As with the simulated point-source, we take this map, smooth and degrade 
it to iVgide values of 256, 128 and 64, and then mask using only the Galactic portion of the 7-year KQ85 mask. The 
resultant maps are plotted in Fig. 11 along with the input. The extra smoothing creates a strip, a few pixels wide, of 
contamination around the Galactic mask, typically at a level of 0.3 — 0.4 //K. Scaling this up by a factor of 50 yields 
contaminants of ~ 20 fiK. Although this is an order of magnitude higher than that created by point-sources, it is 
extremely localized, and does not mimic any of the target signals. We conclude that, as with the point-source mask, 
there is no need to extend the Galactic cut when degrading. Our adaptive-resolution algorithm should be robust 
to smoothing-induced contamination, a hypothesis that will be further tested at a later point by processing a null 
simulation containing realistic foreground residuals. 



VII. TESTING THE FORMALISM'S APPROXIMATIONS 



The formalism set out in Section |TT] includes two main approximations: that the likelihood need only be integrated 
over ranges of {8q, </>o} (and, indeed, # C rit) corresponding to patches containing candidate sources, and that correlations 
between data inside and outside these patches can be neglected. The error made in each of these approximations 
will depend on the nature of the source templates. In particular, the accuracy of the first approximation will depend 
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1% V-ILC Degraded 1% V-ILC (N side = 256) 




-1.1 — 0.027 (J.K -1.7 ^— ^— 2.6e-05 



FIG. 11. One percent of the difference between the WMAP V-band and ILC maps provides an estimate of the Galactic 
contamination present in WMAP data products. This map is plotted here (top-left), alongside the residuals produced by 
smoothing and degrading it to iV s id e ,deg = 256, 128 and 64. All low-resolution plots are shown with a degraded KQ85 Galaxy- 
only mask applied. 



on how well the candidate regions cover the underlying sources, and how the likelihood falls off as a function of the 
location and size of the feature. It is also important to note that the approximations will most likely improve as a 
function of the signal-to-noise for the source: if the likelihood is not peaked in any region of parameter space, we are 
not justified in choosing to integrate only over particular regions. 

To explore these issues, consider the single-source contribution to the posterior Eq. [lO] Neglecting correlations 



between regions inside and outside the patch (see Eq. 15 ) corresponds to the assumption that 



tC^dJ < 1 , (47) 



12 



where the template has support in a blob b, and the data consist of pixels in region b outside the blob. In Fig 
we plot the inverse covariance between several positions and the rest of the sky (e.g., a set of rows of the inverse 
covariance matrix) in ACDM using the best-fit WMAP 7-year cosmological parameters, keeping only the first 50 
multipole moments. It can be seen that the inverse covariance is only significant within a disk of radius ~ 15° around 
each of the template pixels. Therefore, we need not retain all of the pixels on the sky. In Fig. [l3j we depict the case 
where the likelihood is peaked for templates inside a region well-contained within the blob (shaded disk). The inverse 
covariance will be significant within a disk (dashed circle) of ~ 15° around each pixel where the template is non-zero 
(black dot). Our approximation neglects correlations between the template and the pixels contained within the dashed 
circle, but outside the blob (i.e., in region b). The exponential will clearly yield a decreasing correction to the integral 
as the size of the blob is increased, becoming vanishingly small when the radius of the blob is ~ 15° larger than the 
size of templates near the maximum of the likelihood. For terms of higher order in N s , the correction is slightly more 
complicated, but for blobs separa ted by a distance greater than 15°, the argument is the same. Another assumption 



we have made is encoded in Eq. 16 that the pixels in region b do not contribute to the inverse of the covariance 
inside region b. Further, if the actual source is well-contained within the blob, the likelihood will presumably peak in 
a region well-contained within the blob (which is the assumption behind our first approximation of integrating over 
region b alone). 

The reasoning set out above provides qualitative support for the approximations that make this analysis feasible. 
To determine quantitatively how good these approximations are, we have performed three numerical tests. The first 
two tests are designed to determine whether correlations with pixels outside the blobs can indeed be discounted; the 
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FIG. 12. The inverse of the pixel-space ACDM correlation function between 8 — and all 8,8 = 72° and all 8, and 8 = 144° 
and all 8 (left to right). It can be seen that C _1 is largest in magnitude over a ~ 30° window around the pixel being correlated. 




FIG. 13. The case where the source is well-contained within the blob. There is a clearly peaked likelihood for templates 
contained within the shaded disk. For the pixel denoted by the black dot, the inverse covariance is significant only within the 
dashed circle. Our approximation neglects correlations with pixels in the hatched region. This approximation is worst for pixels 
on the edge of the template. 



third determines the effects of restricting the position integral to within our candidate blobs. 



A. Tests of neglected correlations 

The ideal test of the effect of neglected correlations would be to perform evidence integrals for simulated collisions 
of varying sizes using covariance matrices ranging from patch-sized to full-sky. Unfortunately, memory restrictions 
mean we can only hold the full-sky covariance matrix in memory for HEALPix resolutions smaller than iV s i<ie = 64. At 
this resolution each pixel is ~ 1° across, so only large collisions are faithfully represented. As with the tests of the 
adaptive resolution analysis, we therefore split the test up, performing two tests: one checking the effect of neglected 
correlations on the largest scales, and the other at smaller scales. 

The first test uses the map simulated for the largest-scale degradation tests, containing a 30° bubble collision 
template placed at (0o,<M = (45.0°, 45.0°) with amplitude z = 8x 1CT 5 (as plotted in Fig. [9j). The simulation 
is passed to the candidate detection stage, and the patch corresponding to the bubble collision is singled out. The 
evidence ratio is then calculated twice at -/V s idc.dcg = 64, first using the covariance matrix calculated for the patch, and 
second using the covariance matrix calculated on the full-sky. For clarity, the test is carried out without masking, and 
the integration limits on the template parameters are kept constant for both runs. The difference between evidence 
ratios returned will indicate the scale of any error induced by neglecting correlations at the largest scales. 

The second test reuses the smallest map considered in the degradation tests (again, plotted in Fig. [9]), containing a 
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0patch(°) 


Apix 


log/3 


60.1 


12619 


13.1 ±0.3 


180.0 


49149 


13.3 ±0.3 



TABLE IV. The evidence ratios obtained when a patch covariance matrix is used versus the full-sky covariance matrix. Note 
that the full-sky covariance matrix does not quite cover the entire sky: three pixels are left out. This is a consequence of the 
patch-based nature of the algorithm, and does not affect the conclusions. 



flpatch(°) jVpi x logp 

15.4 14670 13.0 ±0.3 

20.0 24167 13.1 ±0.3 

25.0 37388 13.2 ±0.3 

30.0 53338 13.0 ±0.3 



TABLE V. The evidence ratios obtained when the size of the patch covariance matrix is incrementally increased until all 
correlations are included. 



7° bubble collision template with the same amplitude and position as in test one. As in the first test, the candidate 
detection algorithm is applied, and the feature containing the template is extracted. The evidence is calculated first 
using the standard patch size 15.4° in radius, then using progressively larger patches of sky until the patch is 30° in 
radius. At this point, Fig. 12 implies that the covariance matrix contains all pixels significantly correlated with those 
in the candidate collision region. The difference between evidence ratios will therefore indicate the errors associated 
with neglecting correlations on smaller scales. As with the first test, the integration limits used in each case are the 
same, as is the resolution (iV s ide,deg — 256) at which the calculation is performed, and no mask is used. 

The results of the two tests are presented in Tables IV and |Vj In each case, increasing the patch size does not 
change the evidence value obtained beyond MultiNest precision. This supports that the assumption that correlations 
outside of the patch can be neglected, and indicates that doing so does not add a significant source of systematic 
error, given the scale of the variations induced by the nested sampler. 



B. Test of localization of likelihood peaks 



The third assumption test is designed to assess the claim that the likelihood is peaked in position space, and that 
the evidence values obtained are unchanged when only the peaks are considered; this amounts to changing the limits 
of integration in Eq. [THJ This test uses the same 7°-bubble collision map as in the second test. This time, the patch 
radius is held constant at 30° , but the set of central positions sampled is incrementally increased until the template 
can be centered anywhere within the entire patch. The nested sampler therefore has access to greater portions of 
the collision environs - by the fourth run it can sample central positions placing the templates entirely outside the 
simulated collision region - and can provide an estimate of how much evidence is discarded by restricting the template 
position to lie within the blob. 



The results of the test are reported in Table VI The evidence ratio is stable, indicating that sampling from a 
larger range of central positions does not affect the outcome. This implies that the likelihood is indeed well-localized, 
supporting the assumption that the evidence integration need only be carried out over a restricted range of positions. 
The likelihood is plotted as a function of the two angular coordinates in Fig. 14 and indeed is very strongly peaked 
about the position highlighted by the candidate detection stage. 



VIII. NULL TEST: ANALYSIS OF WMAP END-TO-END SIMULATION 



The WMAP data contain a number of components that cannot be included in the pixel-space covariance matrix; in 
particular, the foregrounds are not known precisely, and so their subtraction leaves behind unknown, highly anisotropic 
residuals. It is therefore important to apply the Bayesian analysis pipeline on a null dataset containing estimated 
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Center Range / pa tch Center Range (°) logp 



10% 1.4 13.0 ±0.3 

25% 3.5 13.0 ±0.3 

50% 7.0 13.0 ± 0.3 

75% 10.5 13.1 ±0.3 

100% 14.0 12.8 ± 0.3 



TABLE VI. The evidence ratios obtained when the range of collision centers sampled is incrementally increased until the 
collision can be placed anywhere within the patch. 




35 40 45 50 55 30 40 50 60 

So (°) 0o (°) 

FIG. 14. The localization-test likelihood as a function of co-latitude 9o (left) and longitude <f>o (right), plotted on a logarithmic 
scale. Overplotted as dashed lines are the integration limits used in each of the five runs testing how the evidence changes as 
the central positions are sampled from larger regions. These limits correspond to sampling central positions from 10%, 25%, 
50%, 75% and 100% of the patch by angular radius. The likelihood is very strongly peaked in both angular coordinates. 



foreground residuals - and any other potential systematic effects - to determine whether they generate false-positive 
results. The full end-to-end simulation of the WMAP experiment used to set the optimal-filter detection thresholds 
is perfectly suited to the task, as it is created from simulated time-ordered data, including foreground contamination, 
and is processed by exactly the same pipeline as the WMAP observations. 

The raw optimal filter analysis of the WMAP 7-year W-band end-to-end simulation (with the KQ85 mask applied) 
generates 19 bubble collision candidates and 10 texture candidates. Any candidates which are heavily masked are 
discarded, as are candidates whose estimated range of sizes has no overlap with the relevant prior on 9 cri t- Finally, 
any candidates which clearly correspond to a single feature are merged, leaving a set of 12 and 4 bubble collision and 
texture candidates, respectively. The sizes and locations of these candidates as estimated by the optimal filters are 
tabulated in Table IVlTl 

The patch evidence ratios obtained by applying the adaptive-resolution Bayesian evidence calculation to these 
candidates are also presented in Table |VII| As expected, no candidate exhibits strong evidence in favor of either the 
bubble collision or texture model: the maximum evidence ratios are e -5 ' 4 and e -4 , respectively. Merging these 
results produces the posteriors on the global model parameters as plotted in Fig. [15] The posteriors for both models 
are clearly peaked at N s — 0, confirming our prior knowledge that there are no bubble collisions or textures in the 
end-to-end simulation of the WMAP W-band data. The null test therefore indicates that we should not expect 
un-modeled foreground residuals and unknown systematics to generate false positives in the WMAP data. Further, 
the maximum evidence ratios obtained provide indicators of the level of response foregrounds and systematics can 
produce: any features found in the data should exceed these values in order to be considered interesting. 
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ID 


6l C rit range (°) 


0o (°) 


0o(°) 


logp 


1 


2-4 


63.1 


142.0 


-9.68 ±0.3 


2 


3-5 


50.2 


221.5 


-10.87 ±0.3 


3 


6-8 


126.4 


214.8 


-8.75 ±0.3 


4 


8-12 


43.8 


152.7 


-8.57 ±0.3 


5 


10-14 


19.1 


326.3 


-9.21 ± 0.3 


6 


14-18 


139.7 


11.3 


-7.43 ±0.3 


7 


20-26 


137.0 


42.7 


-6.74 ±0.3 


8 


30-40 


56.8 


154.0 


-6.01 ±0.3 


9 


35-45 


118.6 


232.4 


-5.83 ±0.3 


10 


60-70 


102.0 


326.6 


-6.43 ±0.3 


11 


80-90 


86.1 


305.9 


-6.64 ±0.3 


12 


85-90 


53.9 


19.7 


-5.43 ±0.3 



TIT 


#crit range (°) 


vo{ ) 


<po( ) 


logp 


1 


2-5 


50.2 


221.5 


-4.41 ± 0.3 


2 


4-9 


45.8 


153.1 


-6.01 ±0.3 


3 


4-15 


143.2 


3.2 


-3.99 ±0.3 


4 


7-50 


38.7 


94.8 


-6.10 ±0.3 



TABLE VII. The size ranges and locations of the final 12 candidate bubble collisions (left) and 4 candidate textures (right) 
found in the end-to-end simulation of the WMAP experiment, along with their patch evidence ratios. The angular positions 
are related to Galactic longitude, I, and latitude, b, by the transformations I — <f) and b — 90° — 6. 




FIG. 15. The posterior probabilities of the global parameters of the bubble collision (left) and texture (right) models, given the 
end-to-end simulation of the WMAP 7- year W-band. The posterior is plotted as a function of one parameter, N s , for bubble 
collisions, and two parameters, N B and e, for textures. The regions containing 68% and 95% of the posterior probability are 
indicated by the dotted and dashed lines in the bubble collision plot, and as dark and light regions in the texture plot. Both 
posteriors are strongly peaked at zero sources. 



IX. ANALYSIS OF WMAP 7- YEAR DATA 



Our analysis of the WMAP 7-year W-band foreground-reduced temperature map (with the KQ85 mask applied) 
produces a total of 32 bubble collision candidates and 33 texture candidates. The candidates' locations are visually 
inspected, and those which are mostly obscured by the mask are discarded; candidates found to have no overlap with 
the relevant priors on # cr it are likewise cut. Any candidates which are obviously coincident are merged at this point. 
The remaining candidates are then required to also be significant in an optimal filter analysis of the WMAP 7-year 
V-band foreground-reduced temperature map. This simple check requires that each feature is interesting across a 
range of frequencies, indicating that it is not due to foregrounds. This final cut leaves a set of 11 and 12 bubble 
collision and texture candidates respectively. The most probable sizes and locations of these candidates are tabulated 
in Table \VUJ\ and plotted in Fig. [IB} 

Applying the adaptive-resolution evidence calculation to the candidates produces the patch evidence ratios also 
reported in Table ^7TU\ No single candidate is strong enough to claim a detection on its own. However, as demonstrated 
in Refs. [HI [TT], it is possible for a number of weak candidates to favor the addition of relics to ACDM even if their 
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TABLE VIII. The size ranges and locations of the final 11 candidate bubble collisions (left) and 12 candidate textures (right), 
along with their patch evidence ratios. The size ranges tabulated are those derived from the optimal filter analysis; where they 
extend beyond the relevant priors they are truncated before the evidence calculation. The angular positions are related to 
Galactic longitude, I, and latitude, b, by the transformations I = <j> and b = 90° — 9. Note that, although bubble collision candi- 
dates 8 and 9 overlap considerably, their positions and size ranges are sufficiently discrepant to justify processing individually. 




FIG. 16. Top: the estimated sizes and amplitudes of the bubble collision (left) and texture (right) candidates located in 
the WMAP 7-year data by the optimal filters. Bottom: the patches of WMAP 7-year data passed to the Bayesian evidence 
calculation for each of these candidates. The bubble collision plot shows all of the data involved in the evidence calculation 
for each candidate; for clarity, the texture plot only shows the core region of each patch. Note that the plots of the optimal 
filter candidates (top) contain only the estimated contributions due to additional sources, and the temperature ranges therefore 
differ from the plots of the WMAP data (bottom). 
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FIG. 17. The posterior probabilities of the global parameters of the bubble collision (left) and texture (right) models, given 
the WMAP 7-year data. The posterior is plotted as a function of one parameter, N a , for bubble collisions, and two parameters, 
N a and e, for textures. The most probable regions containing 68% and 95% of the posterior probability are indicated by the 
dotted and dashed lines in the bubble collision plot, and as dark and light regions in the texture plot. Both posteriors are 
strongly peaked at TVs = 0. 



individual evidence ratios are less than one: only by combining the results obtained for all candidates can the overall 
predictive power of the underlying model be revealed. The posteriors on the global parameters of the bubble collision 
and texture models, derived by combining the results from the candidates, are plotted in Fig. |17| both posteriors 
are peaked at zero sources. The texture model's dimensionless scale of symmetry breaking is constrained to be 
2.6 x 1CT 5 < e < 1.0 x 1(T 4 (at 95% confidence), which, as the prior is defined only within the range 2.5 x 10 5 < 
e < 1.0 x 10 -4 , indicates that the WMAP data do not provide any interesting constraint on this parameter. 

The WMAP 7-year data do not favor the addition of either bubble collisions or textures to ACDM. As none of 
the candidates exhibits significant evidence for the addition of sources to ACDM, we do not check the candidates for 
foreground residuals. 



X. DISCUSSION 



In Refs [HI [TU] and [TT] , searches for bubble collisions and textures using earlier versions of the Bayesian source 
detection pipeline were published. Each previous analysis shares a number of candidate features in common with the 
current analysis, allowing consistency checks to be carried out between versions of the pipeline. Comparing results 
between versions is non-trivial, and must take into account each change made to the algorithm. In particular: 

1. The prior on the bubble collision size has changed from uniform in the range 2-11.25° to being proportional to 
sin^crit in the range 2 — 90°. Ceteris paribus, this will reduce evidence ratios previously reported for bubble 
collision candidates, particularly those at small scales. 

2. The bubble collision template previously allowed for a discontinuity at the template boundary with amplitude 
z C rit . This parameter is now set to zero due to updated theoretical results [331 [55] , and the bubble collision 
model considered in this analysis is consequently nested within the model considered previously. The effects 
of removing the edge can be determined exactly using the Savage-Dickey Density Ratio [50]: the change in 
evidence will be the ratio of the posterior and prior probabilities of the edge amplitude, evaluated at z cr i t = 
using the results of the previous analysis, i.e., 



A1 , Pr(z crit |d,old) 

A log p = log 



Pr(z cri t,old) 



rlt=0 



As there was little evidence to support the edge parameter in the earlier analysis, the ratio of posterior to prior 
at z cr it = is typically ~ 10, and the new evidence ratios are boosted accordingly. 

The new analysis replaces the WMAP 7-year KQ75 mask used in all prior analyses with the KQ85 mask, 
revealing ~ 8% more of the sky. The change in the fraction of the sky available to the algorithm increases the 
prior volume on observable source positions by ~ 8% as well, and the log evidence ratios hence decrease by a 
similar amount. 
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TABLE IX. The expected and observed differences in evidence ratio found for the four best bubble collision and texture 
candidates between the previous incarnation of the pipeline and the present analysis. Included in the expected change are the 
new form of the (9 C rit prior for the bubble collisions, the removal of the edge from the bubble collision template, and the change 
in the mask. Note that the evidence ratios and IDs of the texture candidates were not previously published. 



4. The candidate detection method has changed from needlets to optimal filters, and the ranges of size and 
position deemed significant therefore also change. The integration limits for the patch evidences can differ by 
small amounts if the new ranges either reveal or truncate regions of non-zero likelihood. 

In addition, since the first bubble collision analysis we have used an increased number of MultiNest live points 
and tighter tolerance and efficiency settings. The current settings were chosen to ensure accurate calculation of the 
evidence; however, tests show that there is no significant difference in our results due to the new settings. 

Tabic |LX| shows the expected and observed changes in evidence ratio obtained for the four best bubble collision and 
texture candidates processed in the previous analyses. In the majority of cases the differences between the observed 
and expected changes in evidence ratios are consistent to MultiNest precision, but there are two bubble collision 
cases where the new results show a significant difference. 

The first case is the Cold Spot, candidate 3. The memory restrictions of the previous analysis [SIHO] required that 



the 6> cr it range sampled be truncated (compare Tables VII in Ref. [5] and VIII in the current analysis). The new 
adaptive-resolution algorithm allows the full range of 9 cri t estimated by the candidate detection stage to be sampled, 
revealing an additional peak in the posterior and boosting the evidence accordingly. The second case is candidate 4. 
The patches of data used to calculate the evidence for this candidate contain a ~ 4° x 3° region masked by the KQ75 
mask but not by the KQ85 mask. The improvement in evidence ratio most likely derives from uncovering these extra 
pixels, which produce an extra hot contribution to an already hot feature. 

The only differences between the current texture analysis and that of Ref. [TT] are the candidate detection algorithm 
and the mask employed. The candidates that contributed most significantly to the posterior in the previous analysis 
are all still detected by the optimal filters, and the changes in evidence observed are consistent within MultiNest 
precision, given the new mask and small differences in integration limits from changing the candidate detection stage. 



Indeed, the posteriors produced by the two analyses (Fig. 2 of Ref. [TT] and the right-hand plot of Fig. 17) arc almost 
indistinguishable by eye. 



XI. CONCLUSIONS 



We have presented a hierarchical Bayesian algorithm for the detection of spatially-localized sources in high-resolution 
CMB datasets. The algorithm uses the posterior over the global parameters describing the population of sources to 
determine whether their presence is warranted by the data and prior theoretical knowledge. To cope with the volume 
of data available, a conservative approximation to the posterior is calculated by selecting the most likely candidate 
sources using optimal filters and assuming that the remaining data do not contribute to the likelihood. Candidates 
are processed at the highest data resolution possible, given their size and the computing power available. The effects 
of the approximations and adaptive-resolution analysis have been quantified using a suite of tests, and are found to 
be comparable to the typical variance in sampling from the un-approximatcd posterior. 

As a demonstration, the pipeline has been applied to search for evidence of bubble collisions and cosmic textures 
in the WMAP 7-year data. This work removes the size restriction imposed by memory constraints on the previous 
bubble collision analysis [9] [TU] , as well as optimising the candidate detection stages of previous bubble collision and 
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texture analyses [9"HTT] . The WMAP data do not favor the addition of either bubble collisions or cosmic textures to the 
ACDM model: even though such sources provide higher-likelihood fits, they are not sufficiently predictive to overcome 
the extra model complexity. In the context of these models our results also place limits on the average numbers of 
bubble collisions and textures per CMB sky, which are constrained to be fewer than 4.0 and 5.2 at 95% confidence, 
respectively. The WMAP data do not place any significant constraint on the dimensionless scale of symmetry breaking 
for textures, e. 

The Planck satellite |51) will soon release temperature data with a factor of 2-3 improvement in resolution and 
~ 10 in pixel noise over WMAP. Further, it will extract essentially all of the information from the temperature 
power spectrum, providing a near-ideal characterization of the dominant source of noise in the analysis. These facts 
strongly motivate performing the texture and bubble collision analyses on the Planck data when they become available. 
In addition, high-quality CMB polarization data are being gathered by experiments such as Planck, ACTPol 52J, 
SPTPol [53] and Spider [54]. Textures do not induce a polarization signal [21], but bubble collisions are expected to 
create characteristic imprints |27| . complementary to those in the temperature data. Extension of the hierarchical 
Bayesian analysis pipeline to process polarization data, either in isolation or by cross-correlating with temperature 
maps, therefore represents a promising avenue for future tests of these models. Further, such an upgraded pipeline 
could be readily applied to other localized signatures in the CMB, such as the Sunyaev-Zel'dovich effect produced by 
clusters of galaxies. 
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